# Current, as of Thursday, May 28, 2020.
# rm(list = ls())
# setwd("~/Dropbox/Base Dropbox/List Experiment Paper (PS)/Conditional Accept/Replication Files (Carmines-Schmidt PS)/2008 CCAP/Pooled Sensitive Items")
library(foreign)
library(car)
library(readstata13)
library(zeligverse)
library(ggplot2)
#-----
# Input the data.
our.data = read.dta13("CarminesSchmidtPS-replication-2008.dta", generate.factors = TRUE)
names(our.data)
#-----
# Subset the data to remove individuals NOT assigned to one of the five treatment conditions.
our.data = subset(our.data, treatment!="")
table(our.data$treatment)

# Create a new variable indicating whether the respondent was assigned to ANY of the sensitive-item conditions. 
our.data$sensitive.item = ifelse(our.data$treatment!="Control", "Sensitive-Item", "Control")
table(our.data$sensitive.item)
table(our.data$sensitive.item, our.data$treatment)
# Dichotomous, numeric version of this variable (1 = sensitive-item conditions; 0 = treatment condition).
our.data$sensitive.item.binary = ifelse(our.data$sensitive.item=="Sensitive-Item", 1, 0) 
table(our.data$sensitive.item.binary)
table(our.data$sensitive.item.binary, our.data$treatment)
#-----
# Model: treatment effects for full sample, bivariate relationship. 
our.model = lm(item_count ~ sensitive.item.binary, data = our.data)
summary(our.model)

# Simulate: difference-in-means for full sample, bivariate relationship. 
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary, model = "normal", 
              data = our.data)
z.out
# Respondents in the BASELINE condition. 
x.low = setx(z.out, 
             sensitive.item.binary = 0)
# Respondents assigned one of the SENSITIVE-ITEM conditions. 
x.high = setx(z.out, 
              sensitive.item.binary = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals representing the treatment effect, 
# estimated by the model. 
LB.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model: treatment effects as function of PARTISANSHIP. 
our.model = lm(item_count ~ sensitive.item.binary + republican + democrat 
               + sensitive.item.binary*republican + sensitive.item.binary*democrat, 
               data = our.data)
summary(our.model)

# Simulate, treatment effects for REPUBLICANS. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for REPUBLICANS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# REPUBLICANS in the BASELINE CONDITION.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 1,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0)
# REPUBLICANS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 1,
              democrat = 0,
              interaction1 = 1,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for REPUBLICANS. 
LB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate, treatment effects for DEMOCRATS. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for DEMOCRATS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# DEMOCRATS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 0,
             democrat = 1,
             interaction1 = 0,
             interaction2 = 0)
# DEMOCRATS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 0,
              democrat = 1,
              interaction1 = 0,
              interaction2 = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for DEMOCRATS. 
LB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate, treatment effects for INDEPENDENTS. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for INDEPENDENTS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# INDEPENDENTS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 0,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0)
# INDEPENDENTS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 0,
              democrat = 0,
              interaction1 = 0,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for INDEPENDENTS. 
LB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model: treatment effects as function of IDEOLOGY.
our.model = lm(item_count ~ sensitive.item.binary + conservative + liberal 
               + sensitive.item.binary*conservative + sensitive.item.binary*liberal, 
               data = our.data)
summary(our.model)

# Simulate: treatment effects for CONSERVATIVES. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$conservative*our.data$sensitive.item.binary
our.data$interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for CONSERVATIVES.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal 
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# CONSERVATIVES in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 1,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0)
# CONSERVATIVES in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 1,
              liberal = 0,
              interaction1 = 1,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for CONSERVATIVES. 
LB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for LIBERALS. 
# Set up the interaction terms. 
interaction1 = our.data$conservative*our.data$sensitive.item.binary
interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for LIBERALS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal 
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# LIBERALS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 0,
             liberal = 1,
             interaction1 = 0,
             interaction2 = 0)
# LIBERALS in the SENSITIVE-ITEM CONDITIONS.
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 0,
              liberal = 1,
              interaction1 = 0,
              interaction2 = 1)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for LIBERALS.
LB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for MODERATES. 
# Set up the interaction terms. 
interaction1 = our.data$conservative*our.data$sensitive.item.binary
interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for MODERATES.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal 
              + interaction1 + interaction2,  
              model = "normal", data = our.data)
z.out
# MODERATES in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 0,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0)
# MODERATES in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 0,
              liberal = 0,
              interaction1 = 0,
              interaction2 = 0)
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for MODERATES. 
LB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#----------
# Point estimates for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(PE.full, PE.lib, PE.mod, PE.con, PE.dem, PE.ind, PE.gop), digits = 2)
# Lower bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(LB.full, LB.lib, LB.mod, LB.con, LB.dem, LB.ind, LB.gop), digits = 2)
# Upper bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(UB.full, UB.lib, UB.mod, UB.con, UB.dem, UB.ind, UB.gop), digits = 2)
#----------
# Model: treatment effects for full sample, adding confounding variables. 
our.model = lm(item_count ~ sensitive.item.binary + pid_7 + ideology +
                 female + white + black + education + church + rr + obama_favor +
                 interest + income, 
               data = our.data)
summary(our.model)

# Simulate: difference-of-means for full sample, controlling for confounding variables.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + pid_7 + ideology +
                female + white + black + education + church + rr + obama_favor +
                interest + income, 
              data = our.data, model = "normal")
z.out
# Respondents in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             pid_7 = mean(our.data$pid_7, na.rm=TRUE),
             ideology = mean(our.data$ideology, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# Respondents in the SENSITIVE-ITEM CONDITIONS.
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              pid_7 = mean(our.data$pid_7, na.rm=TRUE),
              ideology = mean(our.data$ideology, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for FULL SAMPLE.
LB.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.full = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model treatment effects as a function of PARTISANSHIP, controlling for confounding variables.
our.model = lm(item_count ~ sensitive.item.binary + republican + democrat + 
                 sensitive.item.binary*republican + sensitive.item.binary*democrat +
                 ideology + female + white + black + education + church + rr +
                 obama_favor + interest + income, 
               data = our.data)
summary(our.model)

# Simulate: treatment effects for REPUBLICANS, controlling for confounding variables.
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for REPUBLICANS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat + 
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income,  
              model = "normal", data = our.data)
z.out
# REPUBLICANS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 1,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(our.data$ideology, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# REPUBLICANS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 1,
              democrat = 0,
              interaction1 = 1,
              interaction2 = 0,
              ideology = mean(our.data$ideology, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for REPUBLICANS.
LB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.gop = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for DEMOCRATS. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for DEMOCRATS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat + 
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income,  
              model = "normal", data = our.data)
z.out
# DEMOCRATS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 0,
             democrat = 1,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(our.data$ideology, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# DEMOCRATS in the SENSITIVE-ITEM CONDITIONS.
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 0,
              democrat = 1,
              interaction1 = 0,
              interaction2 = 1,
              ideology = mean(our.data$ideology, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for DEMOCRATS.
LB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.dem = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for INDEPENDENTS. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$republican*our.data$sensitive.item.binary
our.data$interaction2 = our.data$democrat*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for INDEPENDENTS.
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + republican + democrat + 
                interaction1 + interaction2 +
                ideology + female + white + black + education + church + rr +
                obama_favor + interest + income,  
              model = "normal", data = our.data)
z.out
# INDEPENDENTS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             republican = 0,
             democrat = 0,
             interaction1 = 0,
             interaction2 = 0,
             ideology = mean(our.data$ideology, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# INDEPENDENTS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              republican = 0,
              democrat = 0,
              interaction1 = 0,
              interaction2 = 0,
              ideology = mean(our.data$ideology, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for INDEPENDENTS.
LB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.ind = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#-----
# Model: treatment effects as function of IDEOLOGY, adding confounding variables.
our.model = lm(item_count ~ sensitive.item.binary + conservative + liberal +
                 sensitive.item.binary*conservative + sensitive.item.binary*liberal +
                 pid_7 + female + white + black + education + church + rr + 
                 obama_favor + interest + income, 
               data = our.data)
summary(our.model)

# Simulate: treatment effects for CONSERVATIVES. 
# Set up the interaction terms. 
our.data$interaction1 = our.data$conservative*our.data$sensitive.item.binary
our.data$interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for CONSERVATIVES. 
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal +
                interaction1 + interaction2 + 
                pid_7 + female + white + black + education + church + rr +
                obama_favor + interest + income,
              model = "normal", data = our.data)
z.out
# CONSERVATIVES in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 1,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(our.data$pid_7, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# CONSERVATIVES in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 1,
              liberal = 0,
              interaction1 = 1,
              interaction2 = 0,
              pid_7 = mean(our.data$pid_7, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for CONSERVATIVES.
LB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.con = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for LIBERALS. 
# Set up the interaction terms. 
interaction1 = our.data$conservative*our.data$sensitive.item.binary
interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for LIBERALS. 
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal +
                interaction1 + interaction2 + 
                pid_7 + female + white + black + education + church + rr +
                obama_favor + interest + income,
              model = "normal", data = our.data)
z.out
# LIBERALS in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 0,
             liberal = 1,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(our.data$pid_7, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# LIBERALS in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 0,
              liberal = 1,
              interaction1 = 0,
              interaction2 = 1,
              pid_7 = mean(our.data$pid_7, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effects for LIBERALS. 
LB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.lib = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))

# Simulate: treatment effects for MODERATES. 
# Set up the interaction terms. 
interaction1 = our.data$conservative*our.data$sensitive.item.binary
interaction2 = our.data$liberal*our.data$sensitive.item.binary
# Use Zelig to model the difference-in-means for MODERATES
set.seed(47408)
z.out = zelig(item_count ~ sensitive.item.binary + conservative + liberal +
                interaction1 + interaction2 + 
                pid_7 + female + white + black + education + church + rr +
                obama_favor + interest + income,
              model = "normal", data = our.data)
z.out
# MODERATES in the BASELINE condition.
x.low = setx(z.out, 
             sensitive.item.binary = 0,
             conservative = 0,
             liberal = 0,
             interaction1 = 0,
             interaction2 = 0,
             pid_7 = mean(our.data$pid_7, na.rm=TRUE),
             female = median(our.data$female, na.rm=TRUE),
             white = median(our.data$white, na.rm=TRUE),
             black = median(our.data$black, na.rm=TRUE),
             education = median(our.data$education, na.rm=TRUE),
             church = median(our.data$church, na.rm=TRUE),
             rr = mean(our.data$rr, na.rm=TRUE),
             obama_favor = median(our.data$obama_favor, na.rm=TRUE),
             interest = median(our.data$interest, na.rm=TRUE),
             income = median(our.data$income, na.rm=TRUE))
# MODERATES in the SENSITIVE-ITEM CONDITIONS. 
x.high = setx(z.out, 
              sensitive.item.binary = 1,
              conservative = 0,
              liberal = 0,
              interaction1 = 0,
              interaction2 = 0,
              pid_7 = mean(our.data$pid_7, na.rm=TRUE),
              female = median(our.data$female, na.rm=TRUE),
              white = median(our.data$white, na.rm=TRUE),
              black = median(our.data$black, na.rm=TRUE),
              education = median(our.data$education, na.rm=TRUE),
              church = median(our.data$church, na.rm=TRUE),
              rr = mean(our.data$rr, na.rm=TRUE),
              obama_favor = median(our.data$obama_favor, na.rm=TRUE),
              interest = median(our.data$interest, na.rm=TRUE),
              income = median(our.data$income, na.rm=TRUE))
# Generate simulated values. 
s.out2 = sim(z.out, x = x.low, x1 = x.high)
# Extract point estimates and confidence intervals, treatment effect for MODERATES. 
LB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.050))
PE.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.500))
UB.mod = as.numeric(quantile(s.out2$get_qi(qi = "fd", xvalue = "x1")[,1], 0.950))
#----------
# Point estimates for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(PE.full, PE.lib, PE.mod, PE.con, PE.dem, PE.ind, PE.gop), digits = 2)
# Lower bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(LB.full, LB.lib, LB.mod, LB.con, LB.dem, LB.ind, LB.gop), digits = 2)
# Upper bounds for liberals, conservatives, moderates, Democrats, Republicans, independents
# round(cbind(UB.full, UB.lib, UB.mod, UB.con, UB.dem, UB.ind, UB.gop), digits = 2)

